function f = s_radius(ro,ri,r,Z)

n1 = (r-ri).^2+(1+Z).^2;
d1 = (r-ro).^2+(1+Z).^2;
n2 = (r-ri).^2+(1-Z).^2;
d2 = (r-ro).^2+(1-Z).^2;

n3 = (r+ri).^2+(1+Z).^2;
d3 = (r+ro).^2+(1+Z).^2;
n4 = (r+ri).^2+(1-Z).^2;
d4 = (r+ro).^2+(1-Z).^2;

alpha1 = (1+Z).*log((n1./d1).*(n3./d3));
alpha2 = (1-Z).*log((n2./d2).*(n4./d4));

d5 = (r+ri).^2+(1-Z).^2;
d6 = (r-ri).^2+(1-Z).^2;

betai1 = 2*(r+ri).*atan(2*(r+ri)./d5);
betai2 = 2*(r-ri).*atan(2*(r-ri)./d6);

d7 = (r+ro).^2+(1-Z).^2;
d8 = (r-ro).^2+(1-Z).^2;

betao1 = 2*(r+ro).*atan(2*(r+ro)./d7);
betao2 = 2*(r-ro).*atan(2*(r-ro)./d8);

f = alpha1+alpha2+betai1+betai2-betao1-betao2;